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Abstract 

I study a dimer model on the square lattice with nearest-neighbor exclusion as the only interac- 
tion. Detailed simulations using tomographic entropic sampling show that as the chemical potential 
is varied, there is a strongly discontinuous phase transition, at which the particle density jumps 
by about 18% of its maximum value, 1/4. The transition is accompanied by the onset of orienta- 
tional order, to an arrangement corresponding to the {1/2,0,1/2} structure identified by Phares 
et al. [Physica B 409, 1096 (2011)] in a dimer model with finite repulsion at fixed density. Using 
finite-size scaling and Binder's cumulant, the expected scaling behavior at a discontinuous transi- 
tion is verified in detail. The discontinuous transition can be understood qualitatively given that 
the model possesses eight equivalent maximum-density configurations, so that its coarse-grained 
description corresponds to that of the q = 8 Potts model. 
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I. INTRODUCTION 



Lattice models of the fluid-solid phase transition have been studied for many years JH4]. 
Athermal models, characterized by purely excluded volume interactions, are of particular 
interest in this context for their simplicity, given that (off-lattice) systems of hard spheres 
and hard disks exhibit fluid-solid transitions |5|. Analyzing athermal lattice models of this 
kind is relevant not only to the equilibrium melting transition, but to modeling diffusion in 
dense fluids, glassy dynamics, disordered solids, and granular materials. 

In the simplest lattice gas models, the fluid-solid transition is continuous. Thus the model 
with nearest-neighbor exclusion on bipartite lattices belongs to the Ising universality class, 
while the corresponding model on the triangular lattice - the hard-hexagon model - has a 
transition in the q = 3 Potts universality class {(J. Lattice gases with extended exclusion 
regions on the square lattice were recently studied by Fernandes et al. Q]. For exclusion of 
up to n-th neighbors, these authors found, for n — 1 to 5, that the fluid-solid transition is 
discontinuous for n = 3 and 5, and continuous in the other cases. What causes the transition 
to be continuous or not is as yet unclear. It is therefore of interest to study further examples, 
with the goal of identifying which features of a model determine the nature of the phase 
transition. 

Considerable attention has also been given to lattice dimer models, in which each molecule 
occupies a pair of adjacent sites. In the two-dimensional case the model represents adsorption 
of diatomic molecules on a crystalline surface, which has motivated extensive study [9K12]. 
(A complete review of this literature is not feasible here; the interested reader will find 
further references in the papers cited.) Forpurely on-site excluded- volume interactions, it is 
known that no phase transition occurs [3, la]. Inclusion of nearest-neighbor or longer-range 
interactions, whether attractive or repulsive, may lead to the appearance of new phases. In 
addition to transitions between fluid phases of high and low density (as in the monomer 
lattice gas), phases with periodic structure may also arise. Phase diagrams of lattice dimer 
systems have been studied via transfer-matrix approaches 0, n]] and Monte Carlo simulation 
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12|. 



In this work I study a lattice gas composed of dimers that occupy neighboring sites on the 
square lattice, with nearest-neighbor exclusion (see Fig. 1). Using Monte Carlo simulations, 
I show that this model exhibits a strongly discontinuous transition. The high-density phase 



possesses orientational order. To obtain precise results on the phase diagram and thermo- 
dynamic properties, I use tomographic entropic sampling . This Monte Carlo simulation 
method has proved to be efficient and reliable in studies of continuous phase transitions in 
spin models and lattice gases HQ . Thus a further motivation for the present study is to 
assess the utility of tomographic sampling for discontinuous phase transitions. The results 
confirm that the method is indeed useful in this context: a single set of simulations enables 
one to calculate the thermodynamic properties (including the pressure and the absolute 
entropy), as continuous functions of the chemical potential, without metastability effects. 
Detailed comparison with theoretical predictions for finite-size scaling at a discontinuous 
phase transition confirm the reliability of the simulation method, and lend support to the 
idea that the nature of the phase transition is in part determined by the symmetry of the 
close-packed state. 

The balance of this paper is organized as follows. In the next section I define the model, 
followed in Sec. Ill by a description of the simulation method. Sec. IV presents the 
simulation results, while Sec. V contains a discussion of these results and possible directions 
for future investigation. A mean-field analysis using the pair approximation is discussed in 
the Appendix. 



II. MODEL 



Consider a lattice gas consisting of dimers that occupy pairs of nearest-neighbor sites on 
the square lattice, and exclude other dimers from occupying the six nearest neighbor sites. 
An example of a maximum density configuration is shown in Fig. 1. Since the interaction 
between molecules is purely via excluded volume, all allowed configurations have the same 
energy and the model is said to be athermal. Let C denote an allowed configuration, i.e., a set 
of N dimer positions that respect the excluded volume condition. In thermal equilibrium, in 
the grand canonical ensemble, on a lattice of M sites, the probability of this configuration is 
P(C) = e^^/H, where Af) is the grand partition function, and fi denotes the chemical 
potential divided by ksT. (From here on I refer to /j, simply as the chemical potential.) In 
this work I study the model on L x L lattices with periodic boundaries. 

The grand canonical partition function is given by 



Efa,L) = J2 zN V(N,L) (1) 

N=0 

where z = e M and fl(JV, L) is the number of distinct configurations of N molecules on a 
square lattice of L x L sites with periodic boundaries; the maximum number of molecules is 
N m ax — L 2 /4. The pressure is obtained directly from ir = p/ksT = InS/L 2 , while the mean 
number of molecules per site is 

p=^=4e^^ l ) ( 2 ) 

" N 

An analogous formula furnishes (N 2 ), from which we can calculate var(iV)/L 2 , which is 
proportional to the compressibility. One may further define an orientational order parameter, 

m = 52(\ N h ~ N v \) N z N Q(N, L), (3) 

" N 

where (\Nh — N v \)n is the "microcanonical" (fixed- N) average of the absolute difference 
between the numbers of horizontally and vertically oriented molecules. The second and 
fourth moments of this quantity are also studied. 




FIG. 1: A close-packed arrangement of dimers on the square lattice. Since each monomer excludes 
occupancy at its nearest-neighbors, the arrangement forms a hexagonal lattice. 



The simulations reported in Sec. IV indicate that the transition between the low-density, 
unordered phase and the high-density phase is discontinuous. As discussed below, this can 
be understood on the basis of symmetry considerations. In particular, the fact that there 
are eight equivalent close-packed configurations of the kind shown in Fig. 1 suggests that 
this model is equivalent to the 8-state Potts model. Another route to this conclusion arises 
from the observation that the maximal density structure of Fig. 1 has been identified as 
the first of a set of zero-entropy "cus ps" in a dimer model with strongly repulsive, but 
finite, nearest-neighbor interactions {9!. llo|. (These authors denote the structure of Fig. 1 
as {1/2,0, 1/2}.) Varyin g th e temperature at fixed density, p = 1/4, a discontinuous phase 

L The present model corresponds to the zero-temperature limit 



transition is observed 
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of the model studied in 
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lOj (i.e., nearest-neighbor pairs are prohibited). Since the phase 
transition between ordered and disordered phases is discontinuous at finite temperature, it 
should also be so at zero temperature (see Fig. 2). 



ordered 






p=1/4 



FIG. 2: Simplified (schematic) phase diagram of the model with finite repulsion, in the temperature- 
j3p plane. The studies of Refs. [10|] and are realized at density p = 1/4; (3p tends to infinity 
as T — >• along this line. The present study corresponds to T = 0. Since the transition between 
ordered and disordered phases is discontinuous along the line p = 1 /4, it should also be along the 
line T = 0. 



Observe that the orientational order parameter is given by m = L 2 dln'E/dh, where h 
is an external field coupled to the particle orientation; then x — dm/dh is the orientational 



susceptibilty, proportional to var[m]. Note however that m is not the Potts-model order 
parameter, which involves the fractions of dimers in each of the eight close-packed patterns. 
(The latter quantity is not studied in the present work.) Since each close-packed state 
exhibits perfect orientational ordering, it nevertheless seems reasonable to expect m to 
behave similarly to the Potts order parameter. 



III. SIMULATION METHOD 



As is evident from the foregoing discussion, one requires the configuration numbers 
Q(N, L) and the microcanonical averages of \N^ — N v \ in order to evaluate the quantities of 



interest. These are obtained via tomographic entropic sampling [13[], a Monte Carlo method 
in which an initial estimate for the Q(N, L) is refined in a series of iterations, as in the well 
known Wang-Landau sampling (WLS) method [lj]]. As in WLS, I use a target probability 
distribution that is uniform on the set of energy values (or, as in the present case, on the 
set of A values.) Different from WLS, the estimates (denoted Cl(N, L)) are updated only at 
the end of an iteration, using 

n n (N,L) = ^±n n ^(N,L) (4) 

where H n (N) is the histogram at iteration n, i.e., the number of visits to configurations with 
exactly A dimers, and H n is the mean number of visits per class. Each iteration consists in 
ten rather long studies (each involving 2 x 10 7 lattice updates) starting from diverse initial 
configurations. (Three studies begin with a filled lattice, with all dimers oriented vertically, 
three with all dimers oriented horizontally, and four starting from an empty lattice.) The 
final estimates fl are obtained using 5-7 iterations of this kind. The estimates of (\Nh—N v \)^, 
((iV/i — N v ) 2 )n, and ((A/j — A,,) 4 ) ^ are calculated at the final iteration. Further details on 



the method may be found in 13]. 



I study system sizes L — 8, 16, 32, 64, and 128. For the smallest size, I begin with a 
uniform distribution, Cl (N) = 1, VA; after the first iteration the estimates, Qi (A), are 
quite close to their final values. For L — 16, the initial estimate Cl (N) is obtained by 
interpolating the entropy den sity s(p) = \nQ(pL 2 /2,L)/L 2 , using the final estimates in the 



L = 8 study, as explained in 
an analogous manner. 



13(| . The initial estimates for larger systems are generated in 



As noted above, the target probability distribution is uniform on the set of N values, 0, 
1, 2,...,N max . This means that the probability of configuration C, having N dimers, satisfies 
P(C) oc 1/Q(N,L). Using a Metropolis-like procedure, the probability P[C — >■ C] for a 
transition from configuration C to a trial configuration, C, should satisfy, 

P[C -+ C] = n(N,L) 

P[C'^C] Q(N',LY {> 
where N and N' are the numbers of dimers in the current and trial configurations, re- 
spectively. In this way, the transition probabilities satisfy detailed balance with respect 
to the target distribution. Since the Q(N, L) are unknown, the transition probabilities are 
evaluated using the current estimates, Cl. 

The simulation algorithm uses an insertion/removal dynamics. In a removal step, one of 
the N existing dimers is selected at random and deleted. In an insertion step, one of the 
B lattice bonds currently available for insertion is selected and a dimer placed there. (A 
list of available bonds is maintained.) The probabilities for insertion and removal are as 
follows. Let C be a configuration with N dimers and B available bonds. Starting from this 
configuration, insertion is chosen with probability 

p . = i BQ(N, L) 

while removal is chosen with probability 1/2 (independent of the configuration). Unlike 
insertion, however, removal is not always accepted. Suppose the insertion contemplated 
above results in a configuration C . To satisfy the detailed balance relation, Eq. (j3j), a 
removal step yielding configuration C should be accepted with probability 

(N + l)n(N + l,L) 
Pacc {N+l)Q{N+l,L) + BSl{N,L) {> 

This is because the specific transition C — > C then occurs with probability pi ns /B while the 

probability of the inverse transition is p acc /[2(N + 1)]. In general, the acceptance probability 

for removal, starting from a configuration with N dimers, is 

n (MB')- NCl(N, L) 

Pacc{ ' ' NQ(N)+B'Q(N-1,L) {> 
where B' is the number of available bonds in the target configuration. Maintaining of a list 



of available bonds and calculating the change in the number of such bonds before accepting 
a removal involves some computational overhead, but this is more than compensated by the 
improved efficiency, as compared with a "blind" insertion procedure. To choose the next 
transition, a random number z, uniformly distributed on [0, 1) is generated, and if z < Pi ns 
insertion is performed, while if z > 1/2 removal is attempted. In case pi ns < z < 1/2, no 
transition is performed, and the current configuration is counted once again in the histogram 
and the microcanonical averages. 

For each system size, we have the exact values Q(N max ,L) = 8, Q(0,L) = 1, Q(1,L) = 
2L 2 , and Q(2,L) = L 2 (2L 2 — 23). I use the latter value to fix the absolute entropy, via 



^(iV,L) = L 2 (2L 2 -23)^^. (9) 
This is required since the simulations only furnish ratios between configuration numbers. 



IV. RESULTS 

I turn now to results for thermodynamic properties. Figure 3 shows the particle density 
p as a function of chemical potential p, while Fig. 4 shows the pressure versus density, calcu- 
lated directly from the grand canonical partition function. The hallmarks of a discontinuous 
phase transition are immediately evident: a sudden change in density with varying chemical 
potential, associated with a nearly flat region in the isotherm n(p). The orientational order 
parameter m also appears to develop a discontinuity with increasing system size, as shown 
in Fig. 5. From this graph it appears that, in the thermodynamic limit, the low-density 
phase exhibits no orientational order, while in the high-density phase nearly all molecules 
have the same orientation. A striking aspect of these graphs is the absence of metastability 
or hysteresis effects commonly observed in conventional Monte Carlo simulations of a system 
exhibiting a discontinuous phase transition. 

To pinpoint the transition, and elucidate its nature, it is helpful to review some results on 



finite-size scaling at discontinuous phase transitions US ll8H20|. Consider a system that (in 
the thermodynamic limit) suffers a discontinuous phase transition at temperature T*. Let 
T£ c denote the temperature at which the specific heat takes it maximum value, in a system 



of linear size L. The results of Borgs and Kotecky [19[ suggest that at a discontinuous phase 
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FIG. 3: Particle density p versus chemical potential p for system sizes L = 16, 32, 64, and 128, 
in order of increasing steepness. Uncertainties are smaller than the line thickness. Inset: data 
for L = 128 plotted as p = 2(p — ~p)/Ap versus p = p — p* L , compared with the scaling function 
p = tanh&/2. Here ~p is the mean mean of the coexisting densities and Ap their difference; the 
constant b = 347.2. 

transition at which only two phases coexist, the shift, AT£ c = |T* — T£ J, decays ~ L~ 2d in 
d dimensions, for large L. In the case of a discontinuous transition in a g-state Potts model, 



however, the shift decays more slowly 20]: AT[ c ~ L d ; the width of the transition region 
also decreases ~ L~ d . For periodic boundary conditions, the energy and specific heat per site 
follow the scaling forms e ~ (ei + e 2 )/2 - (1/2) (ei - e 2 ) tanh[(ei - e 2 )(/3 - /3*)L d /2 + lng/2] 
and c ~ L d A; B /3 2 /cosh 2 [(ei - e 3 )(/3 - /9*)L d /2 + lng/2]. 

In the present instance, the chemical potential /x is the temperature-like variable, with 
the particle density p playing the role of the energy and K = var(p) that of the specific 
heat. I denote the the chemical potential at which k takes its maximum by pf L K ; that at 
which the variance, x °f the orientational order parameter takes its maximum is denoted by 



/i2 x - Based on the results of 18|, |20|, one expects the shifts, A/i£ = \pi* — p* L \, in chemical 
potential to decrease ~ L~ 2 , and for p and var(p)/L 2 to be governed by scaling functions 
~ tanh/Z and cosh~ 2 7i, respectively, where the scaling variable ~p oc L 2 (p, — p* L ). The inset 
of Fig. 3 confirms the hyperbolic tangent form of the density in the region of the transition. 
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FIG. 4: Pressure ir versus particle density p for system sizes L = 16, 32, 64, and 128, in order of 
increasing flatness in the transition region. Relative uncertainties in ir are of order 10~ 6 . 

As expected, the variances of p and m exhibit sharp maxima in the vicinity of the tran- 
sition; the associated chemical potential values are plotted versus L~ 2 in Fig. 6. Quadratic 
extrapolation (versus L~ 2 ) of the data for L > 16 yields the transition values p* = 4.0076(1) 
(from var[p]) and 4.0077(1) (from var[m]) leading to the estimate p* = 4.00765(15). 

To estimate the densities of the fluid and solid phases at coexistence, I analyze the data 
for p(p) on an interval near p*, but far enough away that the results for the two largest 
system sizes are equal, to avoid effects of finite-size rounding. (The specific ranges are 
3.5 < p < 3.96 and 4.03 < p < 4.5.) Extrapolating these data to p* using a quartic 
polynomial yields coexisting densities pj = 0.1974(1) and p s = 0.2421(1). Applying the 
same procedure to the pressure data, I find ix = 1.00748(2) using the fluid-phase data, 
and 7r = 1.00760(2) from the solid phase; I therefore estimate the pressure at coexistence 
as Hcoex = 1.00754(10). Figure 7 shows a typical configuration in the coexistence region 
(p = 0.22); fluid- and solidlike regions are clearly visible. 

Having determined the coexistence parameters, I return to the scaling properties at the 
transition. Figure 8 confirms that ~K = k/ L d oc 1/ cosh 2 [L d (p — p* L )]. In the case of x = xl L 2 
(Fig. 9), there are significant deviations from the expected scaling, with the maximum 
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FIG. 5: Orientational order parameter m versus chemical potential \x for system sizes L = 16, 32, 
64, and 128, in order of increasing steepness in the transition region. Uncertainties are smaller 
than the line thickness. 



of x apparently growing faster than L 2 . Analysis of the data for L = 16 - 128 yields 
Xmax ~ L 2,18 ^), with some suggestion of a smaller exponent at the largest system size. This 
may be associated with the fact that the jump in m grows with increasing L, as the value 
of m in the fluid phase slowly approaches zero as L grows. It thus appears likely that larger 
systems would be required to observe the scaling of \. In the case of the pressure, the 
deviation from the coexistence value, Ap = p — p CO ex, appears to decrease oc 1/L 2 in the 
coexistence region, as shown in Fig. 10. 

Moment ratios S uch as Binder, ennnnant Q are valuable in analysis of discontinuous 
as well as continuous phase transitions. At a (temperature-driven) discontinuous transition, 



Binder's cumulant, V4 



(E 4 ) /3(E 2 ) 2 exhibits a minimum given by 18|, |20 | 
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ei 



L d 



(10) 



where is a constant and e\ and e 2 are the bulk energy densities of the coexisting phases. 
Identifying, as above, fi as the temperature-like variable, and the density as the variable 
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FIG. 6: Chemical potentials at which var[p] (upper) and var[m] (lower) attain their maxima versus 
1/L 2 . Error bars are smaller than the symbols. 
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FIG. 7: A typical configuration in the coexistence region, p = 0.22, L = 64. 
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FIG. 8: Scaling plot: k = k/L 2 versus \i = L 2 (fi — p£) for (upper to lower) L = 16, 32, 64, and 
128. Relative uncertainties in k are < 4%. Inset: data for L = 128 (points) compared with scaling 
function k = A/ cosh 2 bp, using A = 4.87 x 1CT 4 and b = 0.02182. 



corresponding to the energy density, I define the reduced cumulant as V4 = 1 — (N A ) /3(N 2 ) 2 ; 
Eq. (11 Op is expected to hold if we replace e$ with pj. This is indeed borne out by the 
simulation data: I find (for L = 16 - 128) V A = 0.6531(1) - 0.70(4)/L 2 . Using Eq. (JTUD with 
the estimates for the coexisting densities cited above, I find lim£_ i . 00 V4 = 0.6526(1), quite 
close to the simulation result. The chemical potential values at which V4 takes its minimum 
extrapolate to p = 4.0076(1), consistent with the estimate for the transition point obtained 
previously. 

The results discussed above support, in considerable detail, the conclusion that the model 
exhibits a discontinuous phase transition. Equivalence to the q = 8 Potts model is supported 
by the following observation. The amplitude of the leading correction to the finite size shift, 
i.e., the amplitude A in the relation AT£ c ~ AL~ d , is given by[20| A = (lng)/(ei — e<i). 
Identifying, as before, temperature with chemical potential, and energy density with p, 
we may write, for the effective number of states, q e ff = exp[^4(p! — p 2 )]. The amplitude is 
estimated from the data for \i* L by extrapolating the local slopes, Ah = (A/i^— A{i L , 2 )/[L~ 2 — 
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FIG. 9: Scaling plot: x = xl L 2 versus ~p = L 2 (fi — fi* L ) for (lower to upper) L = 16, 32, 64, and 
128. Relative uncertainties in x are < 4%. 

(L/2)~ 2 ] to L — > oo. The resulting estimate of A = 48.8(6) corresponds to q e ff = 8.0(3), 
consistent with the value expected. 

Finally, it is of interest to note that a discontinuous phase transition is associated with 
a nonconcave portion of the density of states (or microcanonical entropy per site) s(p) = 
lnf2(pL 2 , L)/L 2 . (I employ units such that ks = 1-) The entropy is in general quite difficult 
to obtain using conventional simulation methods, but is accessible via entropic simulation. 
The density of states is compared with the thermodynamic entropy per site, 

s th = ^ 5>(C) lnp(C) = -L L > N ~ (11) 

C N 

in Fig. 11. At the scale of the main graph the microcanonical and thermodynamic entropies 
are indistinguishable, but a detail of s* = s + fi*p shows that the former has a positive 
curvature in the transition region, while the latter, as expected, is concave. 





FIG. 10: Scaled excess pressure it = L 2 (ir — ir coex ) versus density for (lower to upper, on left) 
L = 16, 32, 64, and 128. 

V. DISCUSSION 

A dimer lattice gas with nearest-neighbor exclusion on the square lattice exhibits a 
strongly discontinuous fluid-solid transition accompanied by orientational ordering. To- 
mographic entropic sampling in the grand canonical ensemble proves to be well suited to the 
study of thermodynamic properties, including the absolute entropy and the pressure. The 
expected scaling properties of the Binder cumulant, compressibility and other properties at 
the transition are confirmed, with the possible exception of the variance of the orientational 
order parameter, for which larger system sizes need to be studied. Further analysis of this 
model (and its counterparts on other two- and three-dimensional lattices), should be of con- 
siderable interest, both as regards equilibrium properties, metastability, particle diffusion, 
and nucleation dynamics, as well as nonequilibrium aspects such as a possible glassy phase 
and behavior under an external drive. Adding an attractive interparticle potential, possi- 
bly with orientational dependence, one should expect a richer phase diagram, including a 
liquid- vapor phase transition, nematic ordering, and structural transitions in the solid phase. 

A detailed theoretical analysis of the model appears not to have been performed, although, 
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FIG. 11: Main graph: microcanonical and thermodynamic entropies versus density for L = 128; 
the two functions are indistinguishable at this scale. Inset: s* = s + fi*p for (lower to upper) the 
microcanonical entropy (L = 64 and 128), and the thermodynamic entropy. 



as mentioned in Sec. II, there is strong reason to expect a discontinuous transition given 



the results on the related model with finite repulsion |10l. The transfer- matrix approach 
of Refs. appears most suitable for this analysis. In the Appendix, I show that a 

relatively simple dynamic pair approximation yields a discontinuous phase transition, with 
coexisting densities that are comparable to those found via simulation. 

I close with a few general observations on this and related models. It is worth recalling 
that athermal lattice dimers with on-site exclusion only (the monomer- dimer model, with 
monomers corresponding to sites not occupied by dimers), do not exhibit a phase transition 
[3, Q] • The present study shows that extending the exclusion to first neighbors, a discontinu- 
ous transition arises, in contrast to the simple (monomer) lattice gas, in which exclusion out 
to third neighbors is required for a discontinuous fluid-solid transition. The simple lattice 
gas with nearest-neighbor exclusion exhibits a continuous transition, in the Ising universality 
class, as can be understood on the grounds of symmetry, since there are two symmetric max- 



imally occupied states. For hard hexagons there are three equivalent maximally occupied 



states, and the transition, which is again continuous, belongs to the q = 3 Potts model class 
6] . For the dimer model studied here, there are eight maximally occupied states, suggesting 
equivalence (under a suitable coarse-graining procedure) to the q = 8 Potts model, which 
exhibits a discontinuous phase transition. This correspondence is in fact supported by the 
simulation results. 

It is interesting to note, in this context, that the hard-core lattice gas with exclusion out 
to third neighbors possesses ten distinct maximally occupied configurations. This model also 
suffers a discontinuous phase transition Q|, as would be expected from its similarity to the 
q = 10 Potts model. It seems natural to conjecture, then, that a lattice gas having exactly q 
equivalent maximum-density states (or minimum-energy states) will exhibit a discontinuous 
phase transition if q is sufficiently large, i.e., q > 4 in two dimensions. The converse to 
this statement does not, however, appear to hold. Thus, the number of maximum- density 
configurations of the lattice gas with exclusion out to fifth neighbors is infinite, but the 
model nevertheless exhibits a discontinuous transition js|. 




FIG. 12: Mapping to a particle model: exclusion zone of a particle (central diamond) on sublattice 
A, corresponding to a dimer occupying a horizontal bond on the original lattice. 



Finally, I note that the dimer model may also be mapped to a particle model, i.e., to 
a lattice in which each bond of the original square lattice corresponds to a site on the 
new lattice, which is also square. Under this mapping, an occupied horizontal bond in the 



original model corresponds to an occupied site in sublattice A of the new lattice, while a 
vertical bond corresponds to a site in sublattice B. The set of 22 excluded sites (in the new 
lattice) includes all sites out to fourth neighbors of the central site, as well as two fifth 
neighbors in opposite corners (see Fig. 12); the other two fifth neighbors are not excluded. 
The exclusion zone of a particle in sublattice B is the same, except that the fifth neighbors 
not excluded by a particle in sublattice A are excluded, and vice-versa. The lattice gas with 
exclusion out to fourth neighbors has a continuous fluid-solid transition, while exclusion out 
to fifth neighbors leads to a discontinuous transition js|. The dimer model corresponds to an 
intermediate case, with different orientations of the exclusion zone associated with different 
sublattices. 
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Appendix 



I consider a two-site (pair approximation) mean-field analysis 2l|. This approach seeks 
the equilibrium state of a dynamics in which dimers are inserted at available bonds at rate 
z, and removed at unit rate. To begin, consider a simplified a simplified description. The 
probability that a dimer can be inserted at a given bond is the probability that this bond 
as well as the six neighboring bonds be empty. In the pair approximation this probability 
is estimated as 

\ 6 

q 



Pfree = q f ~J (12) 

where q denotes the fraction of unoccupied bonds and y the fraction of unoccupied sites. 
Let p = N/2L 2 denote the fraction of bonds occupied by dimers. Noting that each dimer 
blocks seven bonds, we may write q — 1 — 7p. Each dimer occupies two sites, so the fraction 
of vacant sites is y = 1 — Ap. (Note that p < 1/8, so that y > 1/2.) Thus we have the 
equation of motion p = z(l — 7p) 7 / (1 — 4p) 6 — p. 

To study the phase transition, the approximation developed above needs to be refined 



somewhat. First, we introduce bond fractions q- L and pi and site fractions ?/, for each sub- 
lattice, i = 1,...,8. In a maximally occupied configuration all the bonds in one sublattice 
are occupied; thus Pi < 1. (Of course, if p^ — 1 for one sublattice, it should be zero for the 
others.) A convenient labeling scheme for the sublattices is shown in Fig. 13. Now, for a 
vacant bond in sublattice j to be available for insertion, the six neighboring bonds, all of 
which belong to other sublattices, must be vacant as well. 

Let Tj be the number of halj -occupied bonds in sublattice j. Each occupied bond in one 
of the neighboring sublattices (i.e., having a shared vertex) contributes to rj; for example, 

ri = P2 +P4 +P5 + Pe +P7 + Ps- (13) 

Knowing the r 3 - and pj, we have for the empty bond fractions, qj = 1 — pj — Tj, while 
Uj — 1 — pj — rj/2. The fraction of occupied bonds in sublattice 1 follows 



yi H2 ha ys y6 m ys n 
z -g pi (14) 



dpi _ qi q 2 g4 ?5 9 6 <?7 Qs 
dt~ Z yf 

with analogous equations for the other sublattices. The set of eight equations for the oc- 
cupation fractions are integrated at fixed z, using a fourth-order Runge-Kutta scheme, un- 
til a steady state is reached. The presence of order is indicated by a nonzero value of 
(ft = Pmax —Pmin, of the difference between the largest and smallest sublattice coverages. 
Starting from an ordered initial distribution (0o > 0.7), the final state is disordered for 
fi < /i< = 1.8378, but takes a nonzero value (<fr > 0.8), for /i > The density jumps 
from p = 0.179 to 0.207 at the transition. The transition point fi < should interpreted as 
a stability limit of the ordered phase ( 1.6. ; clS db spinodal) and not as the coexistence value. 
Indeed, the transition occurs at higher /i values for smaller initial values of 0; for O = 0.6, 
for example, it happens at fi = 1.9023. The coexistence density and chemical potential are 
not furnished by the present dynamic approximation, since it does not provide an estimate 
of the free energy. The dynamic pair approximation nevertheless captures the discontinuous 
nature of the phase transition. 
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FIG. 13: Sublattice labeling scheme used in the pair approximation. 
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